1. Introduction

Circular RNAs (circRNAs) are produced from precursor mRNA (pre-mRNA) back-splicing of thousands of genes in eukaryotes. Although circRNAs are generally expressed at low levels, recent findings have shed new light on their cell type-specific and tissue-specific expression and on the regulation of their biogenesis. Furthermore, the data indicate that circRNAs shape gene expression by titrating microRNAs, regulating transcription and interfering with splicing, thus effectively expanding the diversity and complexity of eukaryotic transcriptomes.

In addition to the circular RNAs that are formed from thousands of loci in eukaryotes by exon back-splicing circularization, different types of circular RNAs exist, which are produced by distinct mechanisms. First, circular viral RNA genomes are ligated by host cellular enzymes to form 3ʹ,5ʹ-phosphodiester bonds or 2ʹ,5ʹ-phosphodiester bonds. Second, circular RNA intermediates can be generated during rRNA processing and permuted tRNA biogenesis in Archaea and algae. Third, a number of housekeeping non-coding RNAs, including small nucleolar RNAs and the ribozyme RNase P, were also identified in circular forms in Archaea. Finally, various circular RNAs that are produced from spliced introns and exons have been reported. For example, excised, self-spliced group I and group II introns are ligated to form circular RNAs with either 3ʹ,5ʹ-phosphodiester bonds or 2ʹ,5ʹ-phosphodiester bonds. Excised tRNA introns are removed by tRNA splicing enzymes and then undergo 3ʹ,5ʹ-phosphodiester ligation in Archaea and Metazoa. Circular RNAs derived from lariat introns have been reported in human cells and in Xenopus tropicalis.

Fig.1 Circular RNA (circRNA) biogenesis and its regulation. Chen,LL. (2016) Nature Reviews Molecular Cell Biology 17[4]:205.

As shown in Fig.1, a. CircRNA is produced by back-splicing circularization, which is catalysed by the spliceosomal machinery. Although it remains unknown exactly how the spliceosome is involved in back-splicing (dashed box), it is proposed that following the sequential assembly of spliceosomal small nuclear ribonucleoproteins (snRNPs) on the precursor mRNA (pre-mRNA), a downstream 5′ donor site of an exon (magenta bars) is joined to an upstream 3′ acceptor site to produce a circRNA. b. Orientation-opposite complementary sequences (red arrows) in introns flanking circularized exons promote circRNA biogenesis by pairing to form RNA duplexes that juxtapose the splice sites. However, alternative splicing may occur and lead to the formation of circRNA with a retained intron (black line). c. Competition for RNA pairing across introns can lead to the formation of different circRNAs from a single gene. d. Muscleblind (Mbl) protein enhances circRNA production from its own gene by binding to introns bracketing the circularized exon15. White diamonds denote Mbl-binding sites. e. A proposed model for quaking (QKI) involvement in circRNA biogenesis during the epithelial–mesenchymal transition (EMT)18. QKI is upregulated during EMT. Low QKI expression in epithelial cells promotes exon inclusion in the linear mRNA, whereas high QKI expression in mesenchymal cells facilitates circRNA formation by juxtaposing the circularized exons. White ovals denote QKI binding sites in QKI-regulated pre-mRNAs. f. A-to-I RNA editing by adenosine deaminase 1 acting on RNA (ADAR1) can suppresses circRNA formation by disrupting RNA pairing between flanking introns.ss, splice site.

2. Materials and Methods

2.1 Library construction and sequencing

Total RNA was extracted using Trizol reagent (Invitrogen, CA, USA) following the manufacturer's procedure. The total RNA quality and quantity were analysis of Bioanalyzer 2100 and RNA 6000 Nano LabChip Kit (Agilent, CA, USA) with RIN number >7.0. Approximately 10 ug of total RNA was used to remove ribosomal RNA according to the manuscript of the Epicentre Ribo-Zero Gold Kit (Illumina, San Diego, USA). Following purification, the ribo-minus RNA fractions is fragmented into small pieces using divalent cations under elevated temperature. Then the cleaved RNA fragments were reverse-transcribed to create the final cDNA library in accordance with a strand-specific library preparation by dUTP method, the average insert size for the paired-end libraries was 300±50 bp. And then we performed the pair-end 2×150bp sequencing on an illumina Hiseq 4000 platform at LC Sciences following the vendor's recommended protocol.

2.2 Bioinformatics analysis

Firstly, Cutadapt[1] and perl scripts in house were used to remove the reads that contained adaptor contamination, low quality bases and undetermined bases. Then sequence quality was verified using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). We used Bowtie2[2] and Tophat2[3] to map reads to the genome of species. Remaining reads (unmapped reads) were still mapped to genome using tophat-fusion[4]. CIRCExplorer[5-6] was used to de novo assemble the mapped reads to circular RNAs at first; Then, back splicing reads were identified in unmapped reads by tophat-fusion and CIRCExplorer. All samples were generated unique circular RNAs. Circular RNA expressions from different samples or groups were caculated by scripts in house. Only the comparisons with p value<0.05 were regarded as showing differential expression by R package edgeR[7].

Bioinformatics pipeline for circRNA (RiboZero)

Reference:

1. Martin M (2011) Cutadapt removes adapter sequences from high-throughput sequencing reads. Embnet Journal 17.

2. Langmead B, Salzberg SL (2012) Fast gapped-read alignment with Bowtie 2. Nature Methods 9: 357-359.

3. Kim D, et al. (2013) TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biology 14: R36.

4. Kim, D., & Salzberg, S. L. (2011). Tophat-fusion: an algorithm for discovery of novel fusion transcripts.Genome Biology , 12(8), R72

5. Zhang, X, et all. (2016). Diverse alternative back-splicing and alternative splicing landscape of circular RNAs. Genome Research, 26(9), 1277-1287.

6. Zhang, X, et al. (2014). Complementary sequence-mediated exon circularization. Cell, 159(1), 134.

7. Robinson, M. D., et al. (2010). Edger: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics, 26(1): 139-40.

3. Project information

3.1 Sample information

Species name: Human

Latin name: Homo sapiens

Specimens: tissue/whole blood

3.2 Group information

3.3 Database

3.4 Bioinformatics software

4. Results

4.1 Quality control and statistics of Reads

Statistics of reads from 6 libraries was shown below. Results of FastQC, PCA analysis and sample correlation were deposited in summary/4_1_Quality_Control.

document location: summary/1_rawdata/ReadsQC.xlsx

4.2 Mapping to genome

First, we utlized Tophat to align all reads on reference genome. Those unmapped reads including back splicing junctions were filtered by Tophat-fusion. Finally, circRNAs were identified by CIRCExplorer.

Statistics of tophat mapping:

Statistics of tophat-fusion mapping:

document location:

summary/2_reference_mapping/2_mapped_stat.xlsx

summary/2_reference_mapping/3_tophat_fusion_mapped_stat.xlsx

Statistics of mapped region:

Statistics of unmapped region (tophat-fusion):

document location:

summary/2_reference_mapping/4_mapped_region_stat.pdf

summary/2_reference_mapping/5_tophat_fusion_mapped_region_stat.pdf

4.3 Localization of circRNA on genome

EM01_M_chromosome.jpegEM02_M_chromosome.jpegEM03_M_chromosome.jpegKTY01_K_chromosome.jpegKTY02_K_chromosome.jpegKTY03_K_chromosome.jpeg

document location: summary/2_reference_mapping/4_DensityDistribution/*_chromosome.jpeg

4.4 Statistics of circRNA typies

Exonic circRNAs were defined as "circRNA". Intronic circRNAs were defined as "ciRNA". Intergenic circRNAs were defined as "intergenic".

document location: summary/3_1_circRNA_candidate_prediction/2_circRNA_circType.jpg

4.5 Overview of circRNA expression

4.5.1 Statistics of circRNA expression

document location:

summary/3_1_circRNA_candidate_prediction/2_circRNA_expression_Statistics.xlsx

summary/3_1_circRNA_candidate_prediction/4_circRNA_expression_Boxplot.png

4.5.2 FPKM Interval (FI) of circRNA expression

document location: summary/3_1_circRNA_candidate_prediction/3_circRNA_expression_Interval.xlsx

4.5.3 Density of circRNA expression

document location: summary/3_1_circRNA_candidate_prediction/5_circRNA_expression_Density.pdf

4.5.4 CircRNA expression profiling

document location: summary/3_1_circRNA_candidate_prediction/1_circRNA_expression.xlsx

4.5.5 Statistics of back splicing junctions

Those reads in sequencing rawdata with back splicing junctions (BSJ) of circRNA could be identified by tophat-fusion and CIRCExplorer. We calculated BSJ of each circRNA and those circRNAs with amount of BSJ were priority to run qRT-PCR.

document location: summary/3_1_circRNA_candidate_prediction/6_circRNA_backsplice_reads.xlsx

4.6 Differential expression analysis of circRNA

4.6.1 Profiling of differentially expressed circRNA

document location: summary/3_2_circRNA_differential_expression/*VS*/1_*VS*_circRNA_differential_expression.xlsx

4.6.2 Volcano of differentially expressed circRNA

document location: summary/3_2_circRNA_differential_expression/*VS*/2_*VS*_Volcano.pdf

4.6.3 Heatmap of differentially expressed circRNA

Heatmap is a graphical representation of data where the individual z values (repeated) or Log10 (non-repeated) contained in a matrix are represented as color scheme. Calculation formula of z value is shown below. Sample-i means expression level (fpkm value).

document location: summary/3_2_circRNA_differential_expression/*VS*/3_*VS*_Heatmap.pdf

4.6.4 GO enrichment analysis of differentially expressed host gene

Gene ontology (GO) (http://www.geneontology.org) is a major bioinformatics initiative to unify the representation of gene and gene product attributes across all species. More specifically, the project aims to: 1) maintain and develop its controlled vocabulary of gene and gene product attributes; 2) annotate genes and gene products, and assimilate and disseminate annotation data; and 3) provide tools for easy access to all aspects of the data provided by the project, and to enable functional interpretation of experimental data using the GO, for example via enrichment analysis.

GO is part of a larger classification effort, the Open Biomedical Ontologies (OBO).

Although gene nomenclature itself aims to maintain and develop controlled vocabulary of gene and gene products, the Gene Ontology extends the effort by using markup language to make the data (not only of the genes and their products but also of all their attributes) machine readable, and to do so in a way that is unified across all species (whereas gene nomenclature conventions vary by biologic taxon).

Significant GO terms were calculated by Hypergeometric equation as shown below. TB gene number=number of total genes; TS gene number=number of differentially expressed host genes of circRNA in total genes; B gene number=total number of gene in all GO terms; S gene number=number of differentially expressed host genes of circRNA in this GO term. Those GO terms with p value<0.05 were defined as significant GO terms.

Significant GO terms::

document location: summary/3_2_circRNA_differential_expression/*VS*/4_GO_Enrichment/*VS*_GO_enrichment_Gene.xlsx

Barplot of enriched GO terms:

document location: summary/3_2_circRNA_differential_expression/*VS*/4_GO_Enrichment/*VS*_GO_enrichment.png

Scatterplot of enriched GO terms:

Number of differentially expressed host genes of circRNA enriched in GO terms, p value and rich factor are shown in scatterplot. Rich factor=(number of differentially expressed host genes of circRNA in GO term)/(total number of genes in GO term). The larger rich factor is, the higher enrichment is.

document location: summary/3_2_circRNA_differential_expression/*VS*/4_GO_Enrichment/*VS*_GO_enrichment_scatterplot.pdf

4.6.5 KEGG enrichment analysis of differentially expressed host gene

KEGG (Kyoto Encyclopedia of Genes and Genomes) (http://www.kegg.jp/) is a collection of databases dealing with genomes, biological pathways, diseases, drugs, and chemical substances. KEGG is utilized for bioinformatics research and education, including data analysis in genomics, metagenomics, metabolomics and other omics studies, modeling and simulation in systems biology, and translational research in drug development.

Significant KEGG pathways were calculated by Hypergeometric equation as shown below. TB gene number=number of total genes; TS gene number=number of differentially expressed host genes of circRNA in total genes; B gene number=total number of genes in KEGG pathways; S gene number=number of differentially expressed host genes of circRNA in this KEGG pathway. Those KEGG pathways with p value<0.05 were defined as significant KEGG pathways.

Significant KEGG pathways::

document location: summary/3_2_circRNA_differential_expression/*VS*/5_KEGG_Enrichment/*VS*_KEGG_enrichment_Gene.xlsx

Scatterplot of enriched KEGG pathways:

Number of differentially expressed host genes of circRNA enriched in KEGG pathway, p value and rich factor are shown in scatterplot. Rich factor=(number of differentially expressed host genes of circRNA in KEGG pathway)/(total number of genes in KEGG pathway). The larger rich factor is, the higher enrichment is.

document location:

summary/3_2_circRNA_differential_expression/*VS*/5_KEGG_Enrichment/*VS*_KEGG_enrichment_scatterplot.png

KEGG pathway maps:

Graphical map objects

The KEGG pathway map is a moleculalr interaction/reaction network diagram represented in terms of the KEGG Orthology (KO) groups, so that experimental evidence in specific organisms can be generalized to other organisms through genomic information. Each map is manually drawn with in-house software called KegSketch, which generates the KGML+ file. This file is an SVG file containing graphics objects that are associated with KEGG objects (see KEGG object identifiers). Basic graphics objects in the reference KEGG pathway maps are:

boxes - ortholog (KO) groups identified by K numbers and, in metabolic maps, reactions identified by R numbers as well

circles - other molecules, usually chemical compounds identified by C numbers, but including glycans identified by G numbers

lines - reactions identified by R numbers in metabolic maps; ortholog (KO) groups identified by K numbers in global metabolism maps

and in organism specific pathway maps that are computationally generated:

boxes - genes or gene products identified by the combination of the KEGG organism code and gene identifiers

These map objects can be searched in the search box at the top of the KEGG PATHWAY page, in the search box in each pathway map, and by the KEGG Mapper tools.

Convention of map number prefix

Each pathway map is identified by the combination of 2-4 letter code and 5 digit number (see KEGG Identifiers). The prefix has the following meaning:

map - Reference pathway

ko - Reference pathway (KO)

ec - Reference pathway (EC)

rn - Reference pathway (Reaction)

org - Organism-specific pathway map

Only the first reference pathway map is manually drawn; all other maps are computationally generated. For metabolic pathways, each box (or line) in the reference map is linked to the K number (KO identifeir), the EC number, and the R number (reaction identifier). The KO, EC, and reaction maps are linked to only one of them. For all metabolic and non-metabolic maps, K numbers are converted to gene identifiers in each organism to generate organism-specific pathways.

As shown above, "map" pathways are not colored, "ko/ec/rn" pathways are colored blue, and organism-specific pathways are colored green, where coloring indicates that map objects exist and are linked to corresponding entries. For global metabolism maps, "map" pathways are fully colored, so that "ko/ec/rn" pathways and organism-specific pathways are generated by reducing the coloring indicating the absence of corresponding entries.

In this study, all KEGG pathway maps were deposited in summary/3_2_circRNA_differential_expression/*VS*/5_KEGG_Enrichment/Pathway_pics/.

4.7 Interaction with circRNA and miRNA

According to some published research articles, circRNA may function as (1) transriptional regulation in nucleus, (2) binding to RNA-binding proteins (RBPs) and RNAs besides miRNAs to form RNA-protein complexes, (3) containing an internal ribosome entry site produced a protein product, (4) regulating mRNA in the cell through limited base pairing and (5) miRNA sponge. We utilized miRanda and Targetscan to predict the interaction with circRNA and miRNA in animal and targetmimics (scripts in house) to predict the interaction with circRNA and miRNA. In animal, we could consider those high confidence results of interaction with the higher Targetscan score and lower miranda free energy.

document location: summary/3_3_circRNA_miRNA_network/1_predicted_miRNAs_bind_to_circRNAs.xlsx

5. Contact us

2575 West Bellfort Street

Suite 270

Houston, TX

77054 USA

Local (713) 664-7087

Toll Free: 1-888-528-8818

Fax: (713) 664-8181

http://www.lcsciences.com/

---